#include"sm.h"
#include"help.h"
#include<fstream>
#include<iostream>
#include<iomanip>
#include<cmath>//abs,sqrt
#include<algorithm>//min,max

using std::ifstream;
using std::abs;
using std::min;
using std::max;
namespace SM
{
  void System_Tz_SM::setupOrbitals()
  {
    //read sp orbitals from file
    ifstream OrbIn(OrbFile.c_str());
    if(!OrbIn) std::cerr<<"can't open Orbital Input file!\n";
    int totalOrbitals;
    OrbIn>>totalOrbitals;
    Orbitals.resize(totalOrbitals);
    int num,l,jj,tz,min,max;
    DataType e;
    for(int i=0;i<totalOrbitals;i++)
      {
	OrbIn>>num>>min>>max>>l>>jj>>tz>>e;
	Orbitals[num-1]=OrbitalType(l,jj,tz,e,min,max);
      }
  }

  const int System_Tz_SM::MatEledim=1;//initialize of static variable
  void System_Tz_SM::setupTwoBodyMat()
  {
    int totalChannels=TwoBodyChannels.size();
    TwoBodyMat.resize(totalChannels);
    for(int channel=0;channel<totalChannels;channel++)
      {
	TwoBodyMat[channel].clear();
	int num2BInChannel=TwoBodyChannels[channel].size();
	TwoBodyMat[channel].resize( (num2BInChannel*(num2BInChannel+1))/2 );
      }
    ifstream MatIn(MatFile.c_str());
    if(!MatIn) std::cerr<<"can't open Mat Input file!\n";
    int totalLines;
    MatIn>>totalLines;

    int J,a,b,c,d;
    TwoBodyMatEle MatEle(MatEledim);
    /*MatEledim=1,the value is interaction*/
    for(int i=0;i<totalLines;i++)
      {
	MatIn>>a>>b>>c>>d>>J;
	for(int j=0;j<MatEledim;j++)
	  MatIn>>MatEle[j];
	--a;--b;--c;--d;//the indice start from zero here, the input interaction orbital order start from 1.
	set2BmatAt(a,b,c,d,J,MatEle);
      }
  }

  void System_Tz_SM::printOrbitals() const
  {
    for(int i=0;i<Orbitals.size();i++)
      {
	std::cout<<Orbitals[i].l<<"\t"<<Orbitals[i].jj<<"\t"<<Orbitals[i].tz<<"\t"<<Orbitals[i].e<<std::endl;
      }
  }


  void System_Tz_SM::printTwoBodyStates() const
  {
    for(int i=0;i<TwoBodyStates.size();i++)
      {
	std::cout<<TwoBodyStates[i].a<<"\t"<<TwoBodyStates[i].b<<"\t"<<TwoBodyStates[i].Tz<<"\t"<<TwoBodyStates[i].Par<<"\t"<<TwoBodyStates[i].J<<std::endl;
      }
  }



  //****************************************************************


  void System_Tz_SM_M::setupOrbitals()
  {
    vector<OrbitalType> OrbitalsTemp;
    vector<int> OrdersTemp;
    //split different m orbitals from System_Tz_SM
    for(int i=0;i<pSystem->Orbitals.size();i++)
      {
	int l=pSystem->Orbitals[i].l;
	int jj=pSystem->Orbitals[i].jj;
	int tz=pSystem->Orbitals[i].tz;
	DataType e=pSystem->Orbitals[i].e;
	for(int mm=-jj;mm<=jj;mm+=2)
	  {
	    OrbitalsTemp.push_back(OrbitalType(l,jj,mm,tz,e));
	    OrdersTemp.push_back(i);
	  }
      }
    //store the Orbitals so that protons and neutrons are separated, proton first.
    //notice that we use this property when generating configurations in Configs.
    //the orbitals with same j but adjacent m is put together, this is important to simplify the cal. of diagonal element of J^2 between configs, see getJ2 in Configs.
    Orbitals.clear();
    Orders.clear();
    Orbitals.reserve( OrbitalsTemp.size() );
    Orders.reserve( OrdersTemp.size() );
    for(int i=0;i<OrbitalsTemp.size();i++)
      {
	if(OrbitalsTemp[i].tz==-1)
	  {
	    Orbitals.push_back(OrbitalsTemp[i]);
	    Orders.push_back(OrdersTemp[i]);
	  }
      }
    maxProtons=Orbitals.size();
    for(int i=0;i<OrbitalsTemp.size();i++)
      {
	if(OrbitalsTemp[i].tz==1)
	  {
	    Orbitals.push_back(OrbitalsTemp[i]);
	    Orders.push_back(OrdersTemp[i]);
	  }
      }
    maxNeutrons=Orbitals.size()-maxProtons;
  }
  
  void System_Tz_SM_M::setupTwoBodyStates()
  {
    TwoBodyStates.clear();
    for(int i=0;i<Orbitals.size();i++)
      for(int j=i+1;j<Orbitals.size();j++)
	{
	  int Tz=(Orbitals[i].tz + Orbitals[j].tz)/2;
	  int Par=(Orbitals[i].l + Orbitals[j].l)%2;
	  int M=(Orbitals[i].mm + Orbitals[j].mm)/2;;
	  TwoBodyStates.push_back(TwoBodyStateType(i,j,Tz,Par,M));
	}
    setupTwoBodyChannels();//noting that the twobody states are sorted after this calling.
  }

  const int  System_Tz_SM_M::MatEledim=System_Tz_SM_M::VSystem::MatEledim;//initialize of static variable
  void System_Tz_SM_M::setupTwoBodyMat()
  {
    setupIndiceIn2BStates();
    int totalChannels=TwoBodyChannels.size();
    TwoBodyMat.resize(totalChannels);
    for(int channel=0;channel<totalChannels;channel++)
      {
	int num2BInChannel=TwoBodyChannels[channel].size();
	TwoBodyMat[channel].clear();
	TwoBodyMat[channel].resize( (num2BInChannel*(num2BInChannel+1))/2 );
	for(int bra=0;bra<num2BInChannel;bra++)
	  for(int ket=bra;ket<num2BInChannel;ket++)
	    {
	      int ab=TwoBodyChannels[channel][bra];
	      int cd=TwoBodyChannels[channel][ket];
	      int a=TwoBodyStates[ab].a;
	      int b=TwoBodyStates[ab].b;
	      int c=TwoBodyStates[cd].a;
	      int d=TwoBodyStates[cd].b;
	      TwoBodyMatEle MatEle(MatEledim);
	      get2BmatOnNeed(a,b,c,d,MatEle);
	      set2Bmat(channel,bra,ket,MatEle);
	    }
      }
  }

  void System_Tz_SM_M::setupIndiceIn2BStates()
  {
    int totalOrbitals=Orbitals.size();
    IndiceIn2BStates.resize((totalOrbitals*(totalOrbitals-1))/2);
    for(int i=0;i<TwoBodyStates.size();i++)
      {
	int a=TwoBodyStates[i].a;
	int b=TwoBodyStates[i].b;
	//a<b
	IndiceIn2BStates[mapab(a,b)]=i;
      }
  }

  void System_Tz_SM_M::get2BmatAt(int a,int b,int c,int d,TwoBodyMatEle & MatEle) const
  {
    MatEle.setZero();
    if(a==b || c==d) return;//pauli principle
    int phase=1;
    //make sure a<b,c<d,(ab)<=(cd)
    if(a>b)
      {
	swap(a,b);
	phase*=-1;
      }
    if(c>d)
      {
	swap(c,d);
	phase*=-1;
      }
    int Mab=(Orbitals[a].mm + Orbitals[b].mm)/2;
    int Mcd=(Orbitals[c].mm + Orbitals[d].mm)/2;

    if( Mab!=Mcd ) return;
    int ab=IndiceIn2BStates[mapab(a,b)];
    int cd=IndiceIn2BStates[mapab(c,d)];
    if(ab==-1||cd==-1) return;
    int channelab=ChannelIndice[ab];
    int channelcd=ChannelIndice[cd];
    //channelab should = channelcd
    if(channelab!=channelcd) return;
    int bra=TwoBodyStateIndice[ab];
    int ket=TwoBodyStateIndice[cd];
    if(bra>ket)
      swap(bra,ket);
    get2Bmat(channelab,bra,ket,MatEle);
    MatEle*=phase;
  }

  void System_Tz_SM_M::get2BmatOnNeed(int a,int b,int c,int d,TwoBodyMatEle & MatEle) const
  {
    MatEle.resize(MatEledim);
    MatEle.setZero();

    if(a==b||c==d) return;//pauli
    int jja=Orbitals[a].jj;
    int mma=Orbitals[a].mm;
    int jjb=Orbitals[b].jj;
    int mmb=Orbitals[b].mm;
    int Tzab=(Orbitals[a].tz+Orbitals[b].tz)/2;;
    int MMab=mma+mmb;
    
    int jjc=Orbitals[c].jj;
    int mmc=Orbitals[c].mm;    
    int jjd=Orbitals[d].jj;
    int mmd=Orbitals[d].mm;    
    int Tzcd=(Orbitals[c].tz+Orbitals[d].tz)/2;
    int MMcd=mmc+mmd;

    int Jmin=max(abs(jja-jjb)/2,abs(jjc-jjd)/2);
    int Jmax=min((jja+jjb)/2,(jjc+jjd)/2);

    if(Tzab!=Tzcd) return;
    if(MMab!=MMcd) return;
    if((Orbitals[a].l+Orbitals[b].l+Orbitals[c].l+Orbitals[d].l)%2) return;//parity check
    if(Jmin>Jmax) return;

    int MM=MMab;
    int ordera=Orders[a];
    int orderb=Orders[b];
    int orderc=Orders[c];
    int orderd=Orders[d];
    VSystem::TwoBodyMatEle Temp;
    for(int J=Jmin;J<=Jmax;J++)
      {
	pSystem->get2BmatAt(ordera,orderb,orderc,orderd,J,Temp);
	if(Temp.empty()) continue;
	if(ordera==orderb) Temp*=sqrt(2.);
	if(orderc==orderd) Temp*=sqrt(2.);
	MatEle+=Temp*cg(jja,jjb,J*2,mma,mmb,MM)*cg(jjc,jjd,J*2,mmc,mmd,MM);
      }
  }


  double System_Tz_SM_M::jplus(int a,int b) const
  {
    if(Orders[a]!=Orders[b])
      return 0;
    if( (Orbitals[a].mm-Orbitals[b].mm)/2 ==1 )
      {
	int jj=Orbitals[b].jj;
	int mm=Orbitals[b].mm;
	double j=0.5*jj;
	double m=0.5*mm;
	return sqrt( j*(j+1)-m*(m+1) );
      }
    else
      return 0;
  }

  double System_Tz_SM_M::jminus(int a,int b) const
  {
    if(Orders[a]!=Orders[b])
      return 0;
    if( (Orbitals[b].mm-Orbitals[a].mm)/2 == 1 )
      {
	int jj=Orbitals[b].jj;
	int mm=Orbitals[b].mm;
	double j=0.5*jj;
	double m=0.5*mm;
	return sqrt( j*(j+1)-m*(m-1) );
      }
    else
      return 0;
  }
  
  void System_Tz_SM_M::printOrbitals() const
  {
    for(int i=0;i<Orbitals.size();i++)
      {
	std::cout<<i<<"\t"<<Orbitals[i].l<<"\t"<<Orbitals[i].jj<<"\t"<<Orbitals[i].mm<<"\t"<<Orbitals[i].tz<<"\t"<<Orbitals[i].e<<std::endl;
      }
  }


  void System_Tz_SM_M::printTwoBodyStates() const
  {
    for(int i=0;i<TwoBodyStates.size();i++)
      {
	std::cout<<TwoBodyStates[i].a<<"\t"<<TwoBodyStates[i].b<<"\t"<<TwoBodyStates[i].Tz<<"\t"<<TwoBodyStates[i].Par<<"\t"<<TwoBodyStates[i].M<<std::endl;
      }
  }

  
  //****************************************************************
  void Configs::setup(int MM,int Par)
  {
    configs.clear();
    //generate partitions first;
    genPartitions();
    vector<int> blocks;//every block has same quantum number except m
    vector<int> shift;//bits should be shifted of every block
    int bits=0;
    for(int i=0;i<Orbitals.size();)
      {
	shift.push_back(bits);
	int jj1=Orbitals[i].jj+1;
	blocks.push_back(jj1);
	bits+=jj1;
	i+=jj1;
      }
    vector< vector<Config> > BlockConfigs(blocks.size());//configs of every block
    for(int i=0;i<partitions.size();i++)
      {
	for(int block=0;block<blocks.size();block++)
	  {
	    BlockConfigs[block].clear();
	    //partitions[i][block] particles in  blocks[block] orbitals
	    Config start,last;
	    for(int j=0;j<partitions[i][block];j++)
	      {
		start[j]=1;
		last[blocks[block]-1-j]=1;
	      }
	    Config config(start);
	    while(true)
	      {
		BlockConfigs[block].push_back(config);
		if(config==last)
		  break;
		config=next(config);
	      }
	  }
	//combine configs between blocks
	vector<Config> configs_partition;//configs in a partition
	configs_partition.push_back(Config(0));
	for(int block=0;block<blocks.size();block++)
	  {
	    vector<Config> temp=configs_partition;
	    configs_partition.clear();
	    for(int k=0;k<BlockConfigs[block].size();k++)
	      for(int j=0;j<temp.size();j++)
		{
		  configs_partition.push_back( (BlockConfigs[block][k]<<shift[block])^temp[j] );
		}
	  }
	for(int j=0;j<configs_partition.size();j++)
	  {
	    if(getMM(configs_partition[j])==MM && getPar(configs_partition[j])==Par)
	      configs.push_back(configs_partition[j]);
	  }
      }
    dim=configs.size();
    // for(int i=0;i<configs.size();i++)
    //   std::cout<<i<<"\t"<<configs[i]<<std::endl;
  }


  void Configs::genPartitions()
  {
    const vector<System_Tz_SM_M::VSystem::OrbitalType> & JOrbitals=pSystem->pSystem->Orbitals;
    int totalOrbitals=JOrbitals.size();
    int ProOrbitals=0,NeuOrbitals=0;//num of proton and neutron J scheme orbitals;
    for(int i=0;i<totalOrbitals;i++)
      {
	if(JOrbitals[i].tz==-1) ++ProOrbitals;
	else if(JOrbitals[i].tz==1) ++NeuOrbitals;
      }
    //generate partitions
    vector<int> n_low_limits(NeuOrbitals),n_up_limits(NeuOrbitals);
    vector<int> p_low_limits(ProOrbitals),p_up_limits(ProOrbitals);
    //setup limits
    int in=0,ip=0;
    for(int i=0;i<totalOrbitals;i++)
      {
	if(JOrbitals[i].tz==1)
	  {
	    n_low_limits[in]=JOrbitals[i].low_limit;
	    n_up_limits[in]=JOrbitals[i].up_limit;
	    //set to default value if input is absurd
	    if(n_low_limits[in]<0)  n_low_limits[in]=0;
	    if(n_up_limits[in]<0)  n_up_limits[in]=JOrbitals[i].jj+1;
	    ++in;	    
	  }
	else if(JOrbitals[i].tz==-1)
	  {
	    p_low_limits[ip]=JOrbitals[i].low_limit;
	    p_up_limits[ip]=JOrbitals[i].up_limit;
	    //set to default value if input is absurd
	    if(p_low_limits[ip]<0)  p_low_limits[ip]=0;
	    if(p_up_limits[ip]<0)  p_up_limits[ip]=JOrbitals[i].jj+1;
	    ++ip;	    
	  }
      }
    vector<Partition> Npartitions,Ppartitions;
    part(vN,NeuOrbitals,Npartitions,n_low_limits.data(),n_up_limits.data());
    part(vZ,ProOrbitals,Ppartitions,p_low_limits.data(),p_up_limits.data());
    //combine neutron proton partitions togather
    partitions.clear();
    for(int i=0;i<Npartitions.size();i++)
      for(int j=0;j<Ppartitions.size();j++)
	{
	  Partition temp(totalOrbitals);
	  for(int k=0;k<ProOrbitals;k++)
	    temp[k]=Ppartitions[j][k];
	  for(int k=0;k<NeuOrbitals;k++)
	    temp[k+ProOrbitals]=Npartitions[i][k];
	  partitions.push_back(temp);
	}
  }

  
  void Configs::setup_simple(int MM,int Par)
  {
    configs.clear();
    //generate proton and neutron separately
    bitset<numBits> start_p,start_n;
    bitset<numBits> last_p,last_n;
    for(int i=0;i<vZ;i++)
      {
	start_p[i]=1;
	last_p[maxProtons-1-i]=1;
      }
    for(int i=0;i<vN;i++)
      {
	start_n[i]=1;
	last_n[maxNeutrons-1-i]=1;
      }
    vector< bitset<numBits> > configs_p;
    vector< bitset<numBits> > configs_n;
    bitset<numBits> config(start_p);
    while(true)
      {
	configs_p.push_back(config);
	if(config==last_p)
	  break;
	config=next(config);
      }
    config=start_n;
    while(true)
      {
	configs_n.push_back(config);
	if(config==last_n)
	  break;
	config=next(config);
	
      }
    //combine proton and neutron configs togather and choose those the total angular momentum projection is MM/2
    for(int i=0;i<configs_n.size();i++)
      for(int j=0;j<configs_p.size();j++)
	{
	  bitset<numBits> temp(configs_p[j]);
	  // for(int k=0;k<maxNeutrons;k++)
	  //   temp[maxProtons+k]=configs_n[i][k];
	  temp^=(configs_n[i]<<maxProtons);
	  if(getMM(temp)==MM && getPar(temp)==Par)
	    configs.push_back(temp);
	}
    dim=configs.size();
    for(int i=0;i<configs.size();i++)
      std::cout<<i<<"\t"<<configs[i]<<std::endl;
  }

  int Configs::getMM(const bitset<numBits> & config) const
  {
    int MM=0;
    for(int i=0;i<maxProtons+maxNeutrons;i++)
      {
	if(config.test(i))
	  MM+=Orbitals[i].mm;
      }
    return MM;
  }

  int Configs::getPar(const bitset<numBits> & config) const
  {
    int Par=0;
    for(int i=0;i<maxProtons+maxNeutrons;i++)
      {
	if(config.test(i))
	  Par+=Orbitals[i].l;
      }
    return Par%2;
  }
  DataType Configs::getME(int i,int j) const
  {
    DataType temp(0);
    bitset<numBits> diff=configs[i]^configs[j];
    bitset<numBits> diffket=diff&configs[j];
    bitset<numBits> diffbra=diff&configs[i];
    int diffCount=diff.count();
    vector<int> pos;
    if(diffCount > 4) return temp;//we only have twobody interaction
    else if(diffCount==4)
      {
	findones(diffket,pos,2);
	int c=pos[1];
	int d=pos[0];
	findones(diffbra,pos,2);
	int a=pos[1];
	int b=pos[0];
	temp=get2BV(a,b,c,d);
	int sign=count(configs[i],a,b) + count(configs[j],c,d);
	temp=sign%2? -temp:temp;
      }
    else if(diffCount==2)
      {
	findones(diffket,pos,1);
	int b=pos[0];
	findones(diffbra,pos,1);
	int a=pos[0];
	int sign=count(configs[i],a,0) + count(configs[j],b,0);
	temp=get1B(a,b);
	for(int k=0;k<maxProtons+maxNeutrons;k++)
	  {
	    if(configs[i].test(k))
	      temp+=get2BV(a,k,b,k);
	  }
	temp=sign%2?-temp:temp;
      }
    else if(diffCount==0)
      {
	for(int k=0;k<maxProtons+maxNeutrons;k++)
	  {
	    if(configs[i].test(k))
	      {
		temp+=get1B(k,k);
		for(int l=0;l<maxProtons+maxNeutrons;l++)
		  {
		    if(configs[i].test(l))
		      {
			temp+=0.5*get2BV(k,l,k,l);
		      }
		  }
	      }
	  }
      }
    return temp;
  }


  double Configs::getJ2(int i,int j) const
  {
    double temp(0);
    bitset<numBits> diff=configs[i]^configs[j];
    bitset<numBits> diffket=diff&configs[j];
    bitset<numBits> diffbra=diff&configs[i];
    int diffCount=diff.count();
    vector<int> pos;
    if(diffCount > 4) return temp;//we only have twobody interaction
    else if(diffCount==4)
      {
	findones(diffket,pos,2);
	int c=pos[1];
	int d=pos[0];
	findones(diffbra,pos,2);
	int a=pos[1];
	int b=pos[0];
	temp=pSystem->jminus(b,d)*pSystem->jplus(a,c)+pSystem->jminus(a,c)*pSystem->jplus(b,d);
	temp-=pSystem->jminus(a,d)*pSystem->jplus(b,c)+pSystem->jminus(b,c)*pSystem->jplus(a,d);
	int sign=count(configs[i],a,b) + count(configs[j],c,d);
	temp=sign%2? -temp:temp;
      }
    else if(diffCount==2)
      {
	findones(diffket,pos,1);
	int b=pos[0];
	findones(diffbra,pos,1);
	int a=pos[0];
	int sign=count(configs[i],a,0) + count(configs[j],b,0);
	for(int k=0;k<maxProtons+maxNeutrons;k++)
	  {
	    if(configs[i].test(k))
	      temp-=pSystem->jminus(k,b)*pSystem->jplus(a,k);
	    else
	      temp+=pSystem->jminus(a,k)*pSystem->jplus(k,b);
	    
	  }
	temp=sign%2?-temp:temp;
      }
    else if(diffCount==0)
      {
	//J-J+ matrix element between configs = sum_{hp}(J-)_hp (J+)_ph , p(h) is orbital (not) occupied by this config. because the orbital with same j but with adjacent  m is adjacent, we simplify to use one loop only
	for(int ii=0;ii<maxNeutrons+maxProtons;ii++)
	  {
	    if(configs[i].test(ii))
	      {
		int jj=ii+1;
		if(jj>=maxNeutrons+maxProtons) continue;
		// for(int jj=0;jj<maxNeutrons+maxProtons;jj++)
		//   {
		if(configs[i].test(jj)) continue;
		double j=0.5*Orbitals[ii].jj;
		double m=0.5*Orbitals[ii].mm;
		temp+=j*(j+1)-m*(m+1);
		//		  }
	      }
	  }
	double M=getMM(configs[i])*0.5;
	//Jz^2+Jz
	temp+=M*(M+1);
      }
    return temp;
  }

  //recursion algorithm to part 'num' particles into numorbs parts.
  void Configs::part(int num,int numorbs,vector<Partition> & partitions,int *low_limits,int *up_limits)
  {
    partitions.clear();
    //divide and conquer
    int numorbs1=numorbs/2;
    int numorbs2=numorbs-numorbs1;
    if(numorbs==0) return;
    if(numorbs==1)
      {
	if(num>=low_limits[0] && num<=up_limits[0])
	  partitions.push_back(Partition(1,num));
	return;
      }
    for(int num2=0;num2<=num;num2++)
      {
	int num1=num-num2;
	vector<Partition> temp1,temp2;
	part(num1,numorbs1,temp1,low_limits,up_limits);
	part(num2,numorbs2,temp2,low_limits+numorbs1,up_limits+numorbs1);
	for(auto part1:temp1)
	  for(auto part2:temp2)
	    {
	      Partition combined_part(numorbs);
	      for(int i=0;i<numorbs1;i++)
		combined_part[i]=part1[i];
	      for(int i=0;i<numorbs2;i++)
		combined_part[i+numorbs1]=part2[i];
	      partitions.push_back(combined_part);
	    }
      }
  }
  
  //****************************************************************


  void SMSolver::diag(int MM,int Par)
  {
    configs.setup(MM,Par);
    int dim=configs.dim;
    std::cout<<"total dim:\t"<<dim<<std::endl;
    if(dim==0)
      {
	eigVals.resize(dim);
	eigVecs.resize(dim,dim);
	J2s.resize(dim);
	return;
      }
    Mat mat(dim,dim);
    for(int i=0;i<dim;i++)
      for(int j=i;j<dim;j++)
	{
	  mat(i,j)=configs.getME(i,j);
	  mat(j,i)=mat(i,j);
	}
    ComplexEigenSolver<Mat> es(mat);
    //    SelfAdjointEigenSolver<Mat> es(mat);
    eigVals=es.eigenvalues();
    eigVecs=es.eigenvectors();

    //norm
    for(int i=0;i<dim;i++)
      {
	DataType norm=eigVecs.col(i).transpose()*eigVecs.col(i);
	eigVecs.col(i)/=sqrt(norm);
      }

    
    //    sort
    quickSort(eigVals,eigVecs,0,dim-1);

    for(int i=0;i<dim;i++)
      for(int j=i;j<dim;j++)
	{
	  mat(i,j)=configs.getJ2(i,j);
	  mat(j,i)=mat(i,j);
	}
    J2s=(eigVecs.transpose()*mat*eigVecs).diagonal();
  }


  void SMSolver::diag_lanczos(int num,int MM,int Par,int step)
  {
    configs.setup(MM,Par);
    int dim=configs.dim;
    std::cout<<"total dim:\t"<<dim<<std::endl;
    if(dim==0)
      {
	eigVals.resize(dim);
	eigVecs.resize(dim,dim);
	J2s.resize(dim);
	return;
      }

    if(num>dim) num=dim;

    SpMat mat(dim,dim);
    vector< Triplet<DataType> > nonzeros;//nonzeros triplet(row,coloum,value)
    for(int i=0;i<dim;i++)
      for(int j=i;j<dim;j++)
    	{
    	  DataType val=configs.getME(i,j);
    	  if(abs(val)>1e-10)
    	    {
    	      nonzeros.push_back(Triplet<DataType>(i,j,val));
    	      if(i!=j)//!!duplicated elements will be summed up by setFromTriplets())
    		nonzeros.push_back(Triplet<DataType>(j,i,val));
    	    }
    	}
    mat.setFromTriplets(nonzeros.begin(),nonzeros.end());
    
    LanczosIter<SpMat> lanczositer(mat);
    LanczosEigenSolver<SpMat> es(lanczositer);
    
    int lanczosdim=num+step;
    while(true)
      {
    	es.compute(lanczosdim-1,false);
    	Vec eigvals=es.eigenvalues();
    	es.compute(lanczosdim,false);

    	double sigma=(es.eigenvalues().head(num)-eigvals.head(num)).norm();
    	//	std::cout<<lanczosdim<<"\t"<<sigma<<std::endl;
    	if(sigma<1e-8)
    	  break;
    	lanczosdim+=step;
      }
          
    
    es.compute(lanczosdim,true);
    eigVals=es.eigenvalues().head(num);
    eigVecs=es.eigenvectors().leftCols(num);
    

    SpMat matJ2(dim,dim);
    nonzeros.clear();//nonzeros triplet(row,coloum,value)
    for(int i=0;i<dim;i++)
      for(int j=i;j<dim;j++)
	{
	  DataType val=configs.getJ2(i,j);
	  if(abs(val)>1e-10)
	    {
	      nonzeros.push_back(Triplet<DataType>(i,j,val));
	      if(i!=j)//!!duplicated elements will be summed up by setFromTriplets())
		nonzeros.push_back(Triplet<DataType>(j,i,val));
	    }
	}
    matJ2.setFromTriplets(nonzeros.begin(),nonzeros.end());
    J2s=(eigVecs.transpose()*matJ2*eigVecs).diagonal();
  }


  void SMSolver::calOccuNums(int orb)
  {
    int totalOrbitals=configs.Orbitals.size();
    OccuNums.resize(totalOrbitals);
    for(int i=0;i<totalOrbitals;i++)
      {
	OccuNums[i]=0;
      }
    
    for(int i=0;i<eigVecs.rows();i++)//for every configuration
      {
	for(int j=0;j<totalOrbitals;j++)
	  {
	    if( configs.configs[i].test(j) ) //if the orbital j is occupied in this configuration
	      OccuNums[j]+=pow(eigVecs(i,orb),2);
	  }
      }
  }
  void SMSolver::printOccuNums()
  {
    using std::cout;
    using std::setw;
    using std::endl;
    cout<<setw(4)<<"l"<<"\t"<<setw(4)<<"jj"<<"\t"<<setw(4)<<"mm"<<"\t"<<setw(4)<<"tz"<<"\t"<<setw(10)<<"OccuNums"<<endl;

    for(int i=0;i<OccuNums.size();i++)
      {
	cout<<setw(4)<<configs.Orbitals[i].l<<"\t"<<setw(4)<<configs.Orbitals[i].jj<<"\t"<<setw(4)<<configs.Orbitals[i].mm<<"\t"<<setw(4)<<configs.Orbitals[i].tz<<"\t"<<setw(10)<<OccuNums[i]<<endl;
      }
  }

  void SMSolver::printStates(int num)
  {
    using std::cout;
    using std::setw;
    using std::endl;
    if(num==-1||num > eigVals.rows()) num=eigVals.rows();//default num=-1, print all
    
    cout<<"Num"<<"\t"<<setw(15)<<"Energy"<<"\t"<<setw(15)<<"J(J+1)"<<endl;    
    for(int i=0;i<num;i++)
      {
	//		cout<<i<<"\t"<<setw(10)<<eigVals(i)<<"\t"<<setw(8)<<int(sqrt(4*J2s(i)+1)-1 + 0.5)/2.<<endl;
	cout<<i<<"\t"<<setw(15)<<eigVals(i)<<"\t"<<setw(15)<<J2s(i)<<endl;
	//	cout<<i<<"\t"<<setw(10)<<eigVals(i)<<"\t"<<setw(8)<<int(J2s(i).real() + 0.5)<<endl;
      }
  }


  
}
